A Markov chain Monte Carlo analysis to constrain dark matter properties with 

directional detection 
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Directional detection is a promising dark matter search strategy. Indeed, weakly interacting 
massive particle (WIMP)-induced recoils would present a direction dependence toward the Cygnus 
constellation, while background-induced recoils exhibit an isotropic distribution in the Galactic rest 
frame. Taking advantage of these characteristic features and even in the presence of a sizeable 
background, it has recently been shown that data from forthcoming directional detectors could 
lead either to a competitive exclusion or to a conclusive discovery, depending on the value of the 
WIMP-nucleon cross section. However, it is possible to further exploit these upcoming data by 
using the strong dependence of the WIMP signal with : the WIMP mass and the local WIMP 
velocity distribution. Using a Markov chain Monte Carlo analysis of recoil events, we show for the 
first time the possibility to constrain the unknown WIMP parameters, both from particle physics 
(mass and cross section) and Galactic halo (velocity dispersion along the three axis), leading to an 
identification of non-baryonic dark matter. 

PACS numbers: 95.35.+d, 14.80.-j 



I. INTRODUCTION 

Directional detection of Galactic dark matter has 
been first proposed by D. N. Spergel flj highlighting 
the fact that even low angular resolution directional 
detectors could be used to show a clear asymmetry in 
the forward/backward distribution of weakly interacting 
massive particle (WIMP) events with respect to the 
direction of the Cygnus constellation. 
Beyond the simple asymmetry feature, it has recently 
been shown that dedicated statistical data analysis of 
forthcoming directional detectors 0-@ could lead either 
to a competitive exclusion Q or to a conclusive discovery 
d, II], depending on the value of the WIMP-nucleon 
cross section. In the latter case, by using a map 
based likelihood analysis and even in the presence of a 
sizeable background, it is possible to show that the main 
incoming direction does correspond to the direction of 
the Cygnus consteUation (£0,60). This is indeed the 
discovery proof of this detection strategy and it has been 
shown that a 10 kg CF4 detector (MIMAC) operated 
during 3 years, would allow for a high significance 
discovery down to a^^ ~ 10"'' pb In this paper, 

we strive to go one step beyond by trying to constrain 
the properties of Galactic dark matter with directional 
detection. 



Indeed, constraining WIMP parameters (mass 



and cross section a„) with upcoming dark matter 
experiments is a main concern of current phenomeno- 
logical studies, using either indirect detection [ll|, 113; 
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direct detection on its own |12h17| . in combination 
with collider data [l9| o r with the measurements of 
halo star kinematics [20| . The quest for a model- 
independent formalism is a difficult task as the signal 
expected in direct detection depends on the properties 
of both the WIMP particle (mass and cross section) 
and the Galactic dark matter halo (three-dimensional 
local WIMP velocity distribution and density). This 
approach is of particular interest in the context of 
competitive upcoming experiments which might be able 
to give positive WIMP detection instead of background 
rejection. M. Drees and C. L. Shan have proposed a 
model-independent reconstruction of the WIMP velocity 
distribution as well as its various moments (mean 
velocity, dispersions, ...), providing the WIMP mass is a 
priori known or deduced from positive signals from 
at least two direct detectors with different target nuclei 
[l6| . The complementary approach is to constrain the 
WIMP properties with the help of a high dimensional 
multivariate analysis and within the framework of a 
general halo model, with a large number of parameters. 
Thus, the main strength of this study, and hence of 
directional detection, is the possibility of constraining 
the properties of both the dark matter particle and the 
dark matter halo with a single experiment. The choice 
of the fitting model must be well motivated e.g. by 
N-body simulations, as it remains as an ansatz. 



Directional detection presents a high identification 
potential thanks to the use of the double-differential 
spectrum d^R/dEj^dflR, also called the directional 
event rate, in a given recoil energy range. Indeed, its 
shape depends both on the WIMP mass and WIMP 
velocity distribution, while the magnitude mainly 
depends on the product of the local WIMP density 
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and the WIMP-nucleon cross section. Within the 
framework of a multivariate recoil event analysis using 
a Markov chain Monte Carlo (MCMC), we show for 
the first time the possibility to constrain, with a single 
directional experiment, the unknown WIMP parameters, 
both from particle physics (m^^ , cr) and Galactic halo 
(velocity dispersion along the three axis), leading to 
an identification of non-baryonic dark matter. It is, of 
course, possible to include external data, e.g. halo star 
kinematics as in [23|, and to relax some astrophysical 
inputs, as po for instance. However, in this work, we 
focus on the contribution of directional detection on its 
own, highlighting the need for future large directional 
detectors. 

The paper is organized as follows. In Sec. II, the dark 
matter halo modeling is introduced while the directional 
detection framework is presented in Sec. III. Then, the 
Markov chain Monte Carlo analysis is detailed in Sec. 
IV, highlighting the performance of such a method in 
the context of directional detection. Sec. V presents the 
results of this 8 parameter analysis for a directional de- 
tector with a sizeable background contamination and in 
the case of a benchmark dark matter model. Departures 
from this input model, by changing the WIMP mass, the 
velocity anisotropy and the background assumptions are 
presented in Sec. VI. 



II. DARK MATTER HALO MODELING 

Direct detection depends crucially on the local WIMP 
velocity distribution [2T| - [23| and it is important to 
investigate the effect of halo modeling on exclusion 
limits and allowed regions. The alternative strategy 
is to build a multivariate analysis, using a halo model 
with a large number of parameters to be constrained by 
the analysis. The isothermal sphere halo model is often 
considered but it is worth going beyond the standard 
assumption especially when considering recent hints in 
favor of triaxiality. 

Indeed, recent results from N-body simulations are in 
favor of triaxial dark matter halos with anisotropic 
velocity distributions and potentially containing sub- 
structures as subhalos (clumps) and dark disk j24| - |27| . 
Moreover, recent observations of Sagittarius stellar 
tidal stream have shown evidence for a triaxial Milky 
Way dark matter halo [l^, with the short axis being 
approximately aligned with the Galactic x axis (toward 
the Galactic center), and the longest with the Galactic y 
axis (in the direction of the solar motion). However, it is 
noteworthy that this result holds true at large radius (60 
kpc) and N-body simulations have shown that there can 
be significant variations of the axis ratios with radius 
[29| . Hitherto, there is no observational evidence of 
triaxiality at solar position. 

The multivariate Gaussian WIMP velocity distribution 



has been first proposed by N. W. Evans et al. [3(1 ■ It 
corresponds to the simplest triaxial generalization of the 
standard isothermal sphere with a density profile p(r) cx 
leading to a smooth WIMP velocity distribution 
without substructure, with a flat rotation curve and in 
dynamical equilibrium. The velocity dispersion tensor 
(Ty, given by the Jeans equations, is symmetric. Thus, 
one can find an orthogonal basis in which the tensor is 
diagonal leading to the following expression of the WIMP 
velocity distribution in the solar system rest frame. 
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(1) 

where the velocity dispersion tensor (Ty — diag[(Ta;, (Ty, az] 
is assumed to be diagonal in the Galactic rest frame and 
vq is the Sun motion with respect to the Galactic rest 
frame. 

The velocity anisotropy /3(r), is then defined as [sTj . 



m = 1 - 



2al 



(2) 



According to N-Body simulations with or without 
baryons [2^ [33 - [3^ . the /? parameter at Rq — 8 kpc 
of the Galactic center, spans the range — 0.4 which is 
in favor of radial anisotropy. 

As stated above, our choice is to develop a high dimen- 
sional multivariate analysis considering a general enough 
halo model, i.e. with a large number of parameters, and 
by constraining all of them with the data analysis of a 
single experiment. The choice of the halo model must 
be carefully done as it remains as an ansatz. Following 
recent results from N-body simulations with baryons 
[25I [35I , the choice of a multivariate Gaussian seems to 
be a reasonable guess, although one could argue that 
deviations are observed in the WIMP velocity distri- 
bution, making it closer to a generalized Gaussian or 
even a double Maxwellian distribution when considering 
the presence of a corotating dark disk. We argue that 
worrying about the exact shape of the WIMP velocity 
distribution seems to be not relevant, in particular when 
taking into account the fact that the resolution of current 
numerical simulations is many orders of magnitude larger 
than the scale of the ultralocal dark matter distribution 
probed by current and future detectors. This is why 
we have chosen a multivariate Gaussian WIMP velocity 
distribution as a fitting model, in a first attempt to 
constrain both the WIMP parameters (m^,a„) and the 
dark matter halo properties using directional detec- 
tion. Effect of nonsmooth halo model with substructures 
and/or streams will be addressed in a forthcoming paper. 

In the following, the input halo model used to gener- 
ate simulated data, is chosen according to two models 
: a standard isotropic halo (/3 = 0) in which case the 
velocity dispersions are linked to the local circular veloc- 
ity vq = 220 km/s as ctx = = CTz = uo/a/2 w 155 
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km/s; and an anisotropic halo (/3 = 0.4), with the fol- 
lowing velocity dispersions {a^ — 200 km/s; = 169 
km/s; (Jy = 140 km/s}. The latter case corresponds to 
the logarithmic ellipsoidal halo model from 's^ with the 
Sun located on the major axis of the halo, the axis ratios 
p and q being equal to 0.9 and 0.8, respectively. This is 
usually taken as an extreme case for the anisotropy, in 
order to avoid instabilities arising when the ratio of any 
of the velocity dispersion is greater than 3. Indeed, as 
discussed in [22, |30] , in order to consider only physically 
relevant model, every velocity dispersions has to satisfy 
the following constraint: ctj/S < ct^ < 3a j. 



into account finite angular resolution required a full cou- 
pling of track reconstruction analysis with this MCMC 
method. As a first step, and as our goal is to show the 
identification potential of directional detection, an ideal 
detector is considered hereafter, i.e. perfect energy and 
angular resolutions. As shown in [7|, the effect of finite 
angular/energy resolution has been shown to be small as 
far as directional exclusion limits are concerned, provid- 
ing the angular resolution is well estimated via detector 
commissioning, e.g. by using a neutron field [36| . 

B. Directional detection 



III. DIRECTIONAL DETECTION 
FRAMEWORK 

A. Detector configuration 

Several dark matter directional detectors are being 
developed and/or operated : MIMAC [1], DRIFT [Ij, 
DM-TPC [| and NEWAGE 16|]. Directional detection 
requires 3D track reconstruction of recoiling nuclei down 
to a few keV with sense recognition. In fine, an ideal 
directional detector should allow one to evaluate the 
double-differential spectrum cPR/dEjidQn in a given re- 
coil energy range [En-^, Efi,^]. The lower bound is due to 
the threshold ionization energy taking into account the 
quenching factor, while the upper bound allows one to 
limit background contamination, as most of the WIMP 
events are concentrated at low-recoil energy. 
In the following, we consider an ideal detector configu- 
ration which could be within reach in a few years. The 
configuration of the MIMAC project is chosen : a 10 kg 
CF4 detector, operated at 50 mbar and allowing 3D re- 
construction of recoiling tracks with sense recognition. 
The chosen recoil energy range is between 5 and 50 keV 
and an exposure f = 30 kg. year is taken into account for 
data simulation. In order to treat realistic cases, we allow 
for a sizeable residual background contamination in the 
data. Indeed, the discrimination of isotropic background 
events from WIMP events has been early recognized as 
the main strength of this detection strategy How- 
ever, as discussed in 0, H, [Tl - fll [H one of the key 
issues for direct detection is the unknown background en- 
ergy distribution. Two extreme cases may be considered 
[l34l^ . [TgI . [TtI : flat or exponentially decreasing with in- 
creasing recoil energy, i.e. with the same feature as the 
WIMP-induced energy spectrum. Within the framework 
of a dedicated statistical data analysis aiming at the iden- 
tification of dark matter, residual background should be 
accounted for and we will show that it does only mildly 
alter the result. 

Energy and angular resolutions are other points to be 
carefully handled. However, it depends on various track 
parameters such as track length, gaz mixture, initial 
track position and direction. A full study of 3D track 
reconstruction is underway [18] and we argue that taking 



The detector velocity in the Galactic rest frame corre- 
sponds to vq, when neglecting the Sun peculiar velocity 
and the Earth orbital velocity about the Sun^52]. We 
consider the value ?70 = 220 km.s~^ along the y axis. 
In such case, the main incoming direction of the WIMP 
signal should be pointing toward {Iq — 90°, 6© = 0°). 
Using the Galactic coordinates {£, b), the WIMP veloc- 
ity is written in the Galactic rest frame as: 

V — u(cos^cos6 X + sin ^ cos 6 y + sin 6 z) (3) 

Following [33], the directional recoil rate is given by. 
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with the WIMP mass, to^ the WIMP-nucleus reduced 
mass, po = 0.3 GeV/c^/cm'^ the local dark matter den- 
sity (see Sec. IV Cl for discussion), the WIMP-nucleus 
elastic scattering cross section, F{Er) the form factor, 
q refers to the recoil direction expressed in the Galac- 
tic coordinates and iimin = ^fn^Ej^ 12771^ is the minimal 
WIMP velocity required to produce a nuclear recoil of en- 
ergy Eft. In the case of an axial coupling and within the 
Born approximation, the expression of the form factor is 
given by [11] : 



F{Er) = 



y/2mNER X R{^X) 



^/2mNER X R{^X) 



(5) 



where R{^X) is the radius of the target nucleus. 
Eventually,^/ (umin, q) is the three-dimensional Radon 
transform [S^l of the WIMP velocity distribution f{v) 
defined as. 



/(t^min,g) = / d^V S{Vmin ~ V.q)f{v) 



(6) 



Geometrically, the Radon transform is the integral of the 
function f{v) on a plane orthogonal to the direction q at 
a distance Vmin from the origin. Using the Fourier slice 
theorem, P. Gondolo found the expression of the Radon 
transform of the multivariate Gaussian to be 32.1 , 
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(7) 
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Together with ([4]), this expression is of particular interest 
in the context of massive MCMC calculations, as it 
allows one to avoid time-consuming evaluation of the 
3D integral of f{v) for each event at each step. It is, 
however, equivalent to the directional recoil rate of [40| . 

In the energy range [Er-^ , Eji^ ] , the expected number of 
WIMP events fig corresponding to a given set of physical 
parameters is given by, 

where ^ is the total exposure. 



IV. THE DIRECTIONAL MARKOV CHAIN 
MONTE CARLO METHOD 

We present a new method, based on a Markov chain 
Monte Carlo analysis, to extract information on dark 
matter from directional data. First, the interest of use of 
MCMC algorithm is outlined. Then, the method is fully 
described in the following section. Discussion on chain 
efficiency will be done in order to prove that MCMC al- 
gorithms are suited for this type of analysis. 



A. Interest of the MCMC algorithm 

As stated above, directional detection offers the pos- 
sibility of using three-dimensional datasets: the recoil 
energy Er and its direction in Galactic coordinates 
(^R, bn). It follows that constraints on dark matter prop- 
erties should be enhanced by the use of directional de- 
tection when compared to directional insensitive detec- 
tion. This is of particular interest when developing a 
high dimensional multivariate analysis aiming at going 
beyond the standard isotropic halo assumption. Indeed, 
the information enclosed in the directional event rate 
{<P R / dE jidQ, b) , i.e. the energy spectrum and the 2D 
shape of the angular distribution, allows one in princi- 
ple to increase the number of degrees of freedom of the 
fitting model. Indeed, adding directional information to 
the energy one allows one to remove degeneracies among 
fitting parameters and hence to deduce consistent con- 
straints. 

In the following, we list the free parameters of our fitting 
model : 

• (ca;, (Jy^ (Jz) the three velocity dispersions of the 
local WIMP velocity distribution, 

• (^o, 6©) referring to the main direction of the re- 
coiling nuclei. It is indeed an unambiguous signa- 
ture of dark matter detection ^] , 

• the WIMP mass. 



• (7„ the WIMP-nucleon cross section directly related 
to (To in the pure proton approximation for the fiu- 
orine target, 

• Rb the background event rate in the considered en- 
ergy range [5, 50] keV. The background events are 
those remaining after the electron/nuclear recoil re- 
jection, based for instance on the length/energy dis- 
crimination (4l| . 

This leads to an eight-parameter analysis of directional 
dark matter dataset allowing us to quantitatively con- 
strain the WIMP properties and the dark matter halo 
profile. Prior ranges are presented in Table HI 

As the number of free parameters is large, grid cal- 
culation of likelihood or functions are not suitable 
due to the exponential growth of the volume of the pa- 
rameter space. Indeed, in order to ensure a scan of all 
the physical parameter space, the regions of interest, i.e. 
the region where the model fits the data, will fill only 
a tiny part of the whole volume. This corresponds to a 
waste of computation time that is avoided using MCMC 
algorithm. Indeed, Markov chains are used in order to 
sample the likelihood (or x^) distribution according to 
Bayesian statistics, enabling the enlargement of the pa- 
rameter space at a minimal computing time cost by fo- 
cusing on the regions of interest. 

In the following we provide a brief description of the 
MCMC, emphasizing its use in the context of directional 
detection. We refer the reader to a more complete de- 
scription, within the framework of cosmic ray physics, in 
1421 and references therein. 



B. Description of the method 

In a general description, an m-dimensional param- 
eter space is described by the following basis 9 = 
5)(2)^ where each element 61*^"^ refers to one 

of the physical parameter of interest. The MCMC al- 
gorithm enables us to sample the conditional posterior 
Probability Density Functions (PDF) of each parameter 
given the data P{d\D), where D refers to the number of 



Parameter 


Prior range 


(GeV/c^) 


(5, 1000) 


logio(cr„ (pb)) 


(-5,-1) 


^0 n 


(-180,-^180) 


be, n 


(-90,-h90) 


<yx,y,z (kin.s"-') 


(5,500) 


(kg-Vear-1) 


(0,50) 



TABLE I: Parameters with their uniform prior ranges used 
for all MCMC analysis. 
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events N, their direction {iR,bii) and their energy Ej^. 
This can be achieved with the use of the Bayes' theorem 
appHed to parameter inference, 



p{e\D) = 



p{D\e) X n(6>) 



(9) 



where P{D) is the data probabihty, also cahed the ev- 
idence, which can be regarded as a normahzation fac- 
tor, n(0) is the prior probabihty indicating the degree of 
behef before observing the data. FinaUy, P{D\6) corre- 
sponds to the hkehhood function written ^{9). In this 
framework, the posterior PDF P{0\D) is the normahzed 
product of the hkehhood function with the priors. Us- 
ing a Bayesian approach, the posterior PDF of each sin- 
gle parameter 9^°'^ is given by the marginalisation of the 
multidimensional P{9\D) distribution over the other pa- 
rameters 




'Steps 



FIG. 1: Correlation function of the 8 parameters from the 
MCMC run applied to our benchmark model (see sec. |V]) 
with the simple multivariate Gaussian proposal function. 



P(6i(")|i3) 



Uf}, V/3e[l,m]\{Q} 



P{9\D) d9^'^\ (10) 



From each one-dimensional PDF, we can estimate the 
expected value of a given parameter and its confidence 
level (CL). The difficulty is then to evaluate the multi- 
dimensional target PDF p{9) = P{9\D). For the above- 
mentioned reasons, instead of using a grid calculation 
algorithm, we developed a MCMC algorithm in order to 
evaluate p{9). This Monte Carlo sampling of the tar- 
get function is done using Markov chains which are a 
sequence of N points in the m-dimensional parameter 
space. 



{0^}^= 



1,...,N 



} 



(11) 



which is constructed according to the Metropolis- 
Hastings algorithm ensuring that the stationary distribu- 
tion of each chain corresponds to the target distribution 
p{9) being sampled. The Metropolis-Hastings algorithm 
is a random walk in the parameter space where each step 
9i^i is derived from the step 9i with the following proce- 
dure : 

• At each step 9i a trial step f?tiiai is generated from 
a proposal distribution q{9tiied\9i). 

• This trial step is accepted or not according to the 
acceptance probability a calculated as follows : 



a = a{9tria.i\0i) = min 1 



P(^trial) g (Atrial I ^») 
piOd g(0l|4ial) 



(12) 



The probability for the trial step to be ac- 
cepted is equal to a. Notice that in the case 
of a symmetric proposal function q we have 
'?(^'triai|^i) = qi0i\9ti-iai), which simplifies the ex- 
pression of eq. (|l3 



• If the trial step is accepted, then 9i^i = 6'triai and 
if not, the chain stagnate at the same point in the 
parameter space leading to 9i+i — 9i. 

Three characteristics of the Markov chains are worth 
being investigated in order to ensure a consistent sam- 
pling of the target function : 

Burn-in length (b): it corresponds to the number of 
steps (or iterations) to be removed from the beginning in 
order to forget the starting point of the random walk. It 
is estimated as the first step reaching the median value 
of the target distribution E[p{9)\ as 



p{9i) > E[p{9)] 



(13) 



Correlation length (I): it is the required minimal length 
between two steps so that they can be considered as un- 
correlated. By construction, each step depends on the 
previous one. Then, in order to get independent steps, 
some subsampling is needed. It corresponds to rejecting 
all steps which are closer than I to each other. The cor- 
relation length Z*^") of each parameter 0^") is estimated 
by computing the autocorrelation function cj"'' where j 
corresponds to the distance between two steps. Indeed, 
Z^"' is defined as the smallest j for which the correla- 
tion function is strictly less than 1/2, i.e cj"'' < 1/2. It 
should be noticed that the limit of 1/2 is arbitrary but 
has been shown to be sufficient in order to consider the 
steps ^l""* and as uncorrelated [l^. Then, the cor- 
relation length of the whole chain I is defined as. 



I 



(14) 



Hence, in order to consider only independent samples 
(steps) ^ind we have subsampled each Markov chain ac- 
cording to the following procedure ^ind = 0i=b+ki with k 



6 



-8686.9 




-8693.4' ' ' — ' ' — ' ' — ^ 10"^' ' ' ' — ' ' — ^ 

10 10' 10 , ^10 , 

Independent samples Independent samples 

FIG. 2: Left panel : mean of the log-likelihood value of each Markov chain E[log(p(&))] as a function of the number of 
independent samples for 10 Markov chains. Right panel : Var[E(p(^))], E[Var(p(^))] and convergence ratio r as a function of 
the number of steps. Both figures come from the MCMC run applied to our benchmark model (see Sec. IVTl . 



being an integer. Figure [T] represents the autocorrelation 
function for the eight parameters from a MCMC analy- 
sis discussed in Sec. |Vl For this chain, we can see that 
the correlation length is equal to 187 due to the WIMP 
mass parameter. Indeed, as explained in the following 
section, the strong correlation between the WIMP mass 
and cross-section will induce larger correlation length. 
Hence, as the correlation length is linked to the stagna- 
tion of the chain, in order to have a smaller value of /, 
the proposal function has to be carefully chosen to ap- 
proximate the target PDF. 

The efficiency of a Markov chain Monte Carlo sampling 
can then be estimated as the fraction of independent sam- 
ples A'ind with respect to the total number of samples N, 
where Nind is given by : 

N -b 

iVind = (15) 

More efficient is a MCMC sampling, since the number of 
rejected samples is lower, leading to a better estimation 
of the target PDF. The quality of the estimation of 
the target PDF is directly affected by the MCMC 
sampling efficiency and hence by the burn-in and 
correlation lengths. Depending on the input values of 
the different parameters used to simulate pseudodata of 
single directional detection experiments, the sampling 
efficiency is between 0.6% and 8%. It mainly depends 
on the correlation lengths which were found to be 
between 6 and 130 using the second proposal function : 
the multivariate Gaussian with covariance matrix (see 
below). The sampling efficiency could be enhanced by 
using other proposal functions not necessarily Gaussian 
like the Binary Space Partitioning first introduced in 
MCMC sampling by A. Putze et al. However, 
as we are running a large number of Markov chains in 
parallel with a low computational time, such efficiencies 
are largely enough to get well sampled PDFs. 



chain convergence: it is a key criteria worth being in- 
vestigated for MCMC sampling as it ensures that the 
target PDF is being sampled by the different chains. In- 
deed, the left panel of Fig.[2]presents the evolution of the 
mean of the log-likelihood value of each Markov chain 
E[log(p(0))] as a function of the number of independent 
samples for 10 Markov chains. From this figure, we can 
appreciate the fact that the mean value of each Markov 
chain is converging to the same value as well as the vari- 
ance, not shown here. Then, in order to quantify the 
convergence, we can form the following ratio : 



Var 




E 


Var(p(^)) 



As seen on the right panel of Fig. EJ Var[E(p(6'))] tends 
to whereas E[Var(p(^^))] tends to a finite number 
leading to r — when increasing the number of inde- 
pendent samples. Then, we can arbitrarily fix a limit rc 
which will correspond to a chain convergence if r < Tc. 
Following fi^, we have chosen rc = 0.2. Then, from 
the Fig. we can see that for this MCMC run with 10 
independent chains, the convergence status is reached 
at the 40*'' independent step (sample). However, even 
if the chain convergence criteria is reached with only 
a few tenth of independent samples, one should have 
longer chains of independent samples to get a very 
precise estimation of the target function. In our case 
we have run between 10 and 100 Markov chains of 10^ 
steps in order to get more than 5 x 10** independent 
samples for each analysis. Indeed, another interest in 
computing several Markov chains in parallel is that we 
can add all the independent samples from every chain 
together to enhance the estimation of the target function. 
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In this paper, we have considered flat prior for each 
parameter {m^ , log^Q (ct„ ) , , , o-j, , , cr^ , i?f,} . In such 
case, the Bayes' theorem is simphfied and the target dis- 
tribution p{9) reduces to the hkehhood function Ji'{9). 
The latter is given by the extended likelihood function 
(see Ref. [11]) dedicated to unbinned data as. 



n 



S[Rn) ■ 



Ms + Mb 



B{R„ 



(17) 



where and fib — Rb x S. are the expected number of 
WIMP events and background events respectively. i?„ 
refers to the energy and direction of each event while the 
functions S and B are the directional event rate of the 
WIMP events and the background events, respectively. 

As previously highlighted, in order to optimize the 
MCMC sampling efficiency, the proposal function must 
be as close as possible to the target PDF. Two different 
successive proposal functions are used : 

• A multivariate Gaussian in the same basis of the 
parameter space with dispersions cr*-"-' taken from a 
fast evaluation of the likelihood function on a grid. 
In such case, we have. 



g(") 
"^trial 



(18) 



where x is a random variable distributed following 
the normal distribution ^(0, 1). 

• A multivariate Gaussian with the covariance matrix 
estimated from the previous run. Then, the next 
step is calculated using: 



Atrial — Oi 



PCx 



(19) 



where C is the eigenvalue of the covariance matrix, 
P is the matrix of the corresponding eigenvectors 
and £ is a vector of m random variables distributed 
following ^(0, 1). 

In both cases, as the proposal function is a Gaussian, 

we are in the case where (7(6'triai|^i) = 9(^i|^triai) which 
simplifies the expression of the acceptance (Eq. [T^]) . 



input model is used to generate simulated data in 
a 10 kg CF4 detector (as proposed by the MIMAC 
collaboration) with a three-year exposition time. These 
data are then analyzed with the directional MCMC 
method fSec. lIVI) . As stated above, the eight parameters 
{■m^,\og^Q{crn),£Q,bQ,ax,ay,az,Rb}, are taken as free 
parameters in the MCMC analysis, with flat priors, thus 
ensuring that the study is model-independent from a 
point of view of both particle physics (WIMP properties) 
and Galactic physics (halo properties). In particular, 
no previous knowledge of the Galactic dark matter 
halo is needed, the goal being to extract the posterior 
PDF of all parameters and check their consistencies 
with the input model. Departure from isotropy as well 
as modification of the various parameters of the input 
model will be studied in Sec. I VII 



FigureOpresents marginalized distributions (diagonal) 
and 2D correlations (off-diagonal) plots of the eight pa- 
rameters of the analysis of simulated data obtained with 
the benchmark input model. The complete result for 
each parameter, as extracted from marginalized distri- 
butions, is summarized in Table HIl where the output pa- 
rameters are characterized by the mean value extracted 
from their ID posterior PDF, while the error bars are 
accounted for a 68% confidence level. However, to fully 
understand correlations between the parameters, the full 
set of 2D correlations is needed. Moreover, to quantify 
those correlations among the eight different parameters, 
the correlation matrix defined as 



^var[6l(")]var[6'('3)] 
is given in Fig. U and will be discussed hereafter. 



(20) 



The result obtained is threefold : the discovery proof 
is given by the reconstruction of the main incoming 
direction (^0,60) (Sec. IV Ap . Then, the three velocity 
dipersions and hence the velocity anisotropy parameter 
(/3) of the dark matter halo are assessed (Sec. IV B[) 
leading to a constraint on the properties of the WIMP 
particle on the (m^, log]^Q(tT„)) plane, within the frame- 
work of our ansatz. In the following, we detail the 
results arising from the analysis of Fig. [31 



V. RESULT FOR A BENCHMARK INPUT 
DARK MATTER MODEL 

For concreteness, we exemplify this directional MCMC 
method by studying the case of a given benchmark input 
model, i.e. the standard isothermal sphere with an 
isotropic velocity distribution with /3 = (see Sec. |TT] 
for more details). A sizeable background contamination 
(10 kg~^year~^) is accounted for, with a flat energy 
spectrum. We consider a 50 GeV/c^ WIMP with a 
WIMP-nucleon axial cross section cr„ = 10"'^ pb. The 



A. Discovery proof 

Following a previous study ^] , we first present the ex- 
traction of the main incoming direction of the events 
(^0,^0), from a pseudodata analysis. This is a blind 
analysis as these two parameters are taken as free param- 
eters of the analysis. It can be concluded from marginal- 
ized distributions of Fig. [3] that the recovered main recoil 
direction is pointing towards the Cygnus constellation 
within 2.5° at 68 % CL, corresponding to a nonambigu- 
ous detection of particles from the Galactic halo which is 
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TABLE II: Comparison of the values of the parameters for the input model and as extracted after the MCMC analysis from 
the marginalized distributions. We quote mean value of the PDF distribution and (68 % CL) error bars. 




FIG. 4: Correlation matrix as defined by Eq. [20]for the eight 
parameters of the MCMC analysis in the case of an isothermal 
halo with a WIMP mass of 50 GeV.c-2 and a WIMP-nucleon 
cross section. The grey scale represents the absolute values of 
p"'''. Signs of correlation can be deduced from Fig. [3] 



in favor of a dark matter positive detection. 
We found the same result as in [8|, which was expected 
as the information on the main incoming direction is en- 
closed mainly in the angular part of the WIMP-induced 
spectrum. The major update is the fact that we prove 
that this result is model-independent as no a priori 
knowledge of neither the Galactic dark matter halo nor 
the WIMP particle is needed. The 2D correlations of 
Fig. [3] (third and fourth columns and lines) and the two 
first columns of the correlation matrix (see Fig. 2]) indi- 
cate that there is no correlation of the main incoming 
direction with the six other parameters. 
In addition, this result holds true for all cases studied 
hereafter where we have considered different input val- 
ues of the WIMP mass or halo model. Indeed, for each 
case, we have checked that the main incoming direction 
is always consistently constrained and reveals no correla- 
tion with other parameters (see Secl VIj) . 
We emphasize conclusions from Q : directional de- 
tection of dark matter is a powerful strategy to clearly 
identify a positive dark matter signal, using the main in- 
coming direction as the discovery proof, even in the case 
of a sizeable background contamination and non standard 
halo model. 



B. Dark matter halo properties 

The originality of this work, in comparison to current 
phenomenological studies, is that the properties of the 
dark matter halo itself are constrained using a single di- 
rectional detection experiment. As shown on Fig. [3l the 
velocity dispersions are strongly and consistently con- 
strained according to the input values. Indeed, from the 
marginalized distributions of the posterior PDF of each 
velocity dispersion, the following constraints can be de- 
duced 

cr^ = IbStlr km.s-i (68% CL), 
ay = 1641^7 km.s-i (68% CL), 
a, ^ km.s-i (68% CL), 

giving, in such case, strong evidence in favor of an 
isotropic dark matter halo. 

However, one could notice, from the full MCMC result 
(Fig. [3]) and from the correlation matrix (Fig. 2]), that 
the three velocity dispersions are quite correlated to 
each other, to the WIMP properties {m y^, log iq {an)) and 
to the background rate. In the following, we propose a 
short discussion to understand the fundamental origin 
of these different correlations. To begin with, the 
positive correlation between each of the three velocity 
dispersions mainly comes from the information on the 
angular distribution. Indeed, in order to reproduce the 
shape of the velocity distribution, which is isotropic 
in this case, the three velocity dispersions have to be 
positively correlated to each other; in this case, we 
found p[ai,aj] w 0.4 for i ^ j e {x,y,z}. However, 
increasing the velocity dispersions leads to an increase 
in the number of expected WIMP events and to 
wider WIMP event angular distribution. The latter 
can be compensated by decreasing the WIMP mass 
as it leads to tighter angular distribution (see 0, Q 
for a detailed discussion) thus implying a negative 
correlation between the WIMP mass and the three 
velocity dispersions with p[m^,ax] — p[m-^,az] ~ —0.55 
and p[m^,ay] ~ —0.75. As the cross section is directly 
proportional to ps the correlations between the log]^g(cr„) 
and the three velocity dispersions are obviously negative 
with p[logiQ{a„),ax] = /9[logio(CT„), cr^] « -0.57 and 
p[logig((T„), (Ty] w —0.70. Correlations between the 
WIMP parameters and ay are stronger than in the case 
of fJa; and (Tz as it is the most related to ps and the total 



10 




-1 -0.5 0.5 



Velocity anisotropy (3 

FIG. 5: Posterior PDF distribution of the /3 parameter, with 
and without correction due to nonflat prior. The prior is 
Monte Carlo estimated. 



width of the angular event distribution. Then, as the 
velocity dispersion along the y axis is more degenerated 
with the other parameters than Ux and az the error bar 
on the estimation of Uy are larger than for the two other 
velocity dispersions (about 2 times larger). The negative 
correlation between the three velocity dispersions and 
the background rate can be easily explained by the 
definition of the extended likelihood function where the 
sum of /is (proportional to the velocity dispersions) 
with /ifc (proportional to the background rate) follows a 
Poisson distribution of mean equal to fXs + fJ-b = -^cvcnt! 
we found in this case: p[Rb^ crj] ~ —0.4 with j e {x,y,z}. 



The evaluation of the velocity anisotropy parameter /3 
allows us to summarize the results from the three velocity 
dispersions. Indeed, the posterior PDF of the /3 param- 
eter can be computed from Eql2] However, having a flat 
prior on the three velocity dispersions implies a nonflat 
(informative) prior on the /3 parameter. Hence, Fig. [5] 
presents the raw PDF of /3 considering flat priors on the 
(Ti's in the black solid line, the induced prior on the /3 
parameter H(/3) is shown as the blue dashed line, and 
the red dotted line corresponds to the corrected PDF of 
/3 with a flat prior. In the following, for each case, we 
will only consider the corrected posterior PDF of the /3 
parameter. From the latter, we can deduce an interest- 
ing constraint /3 = -OmStoAl (68% CL) favoring an 
isotropic dark matter halo. This is a proof that within 
the framework of the multivariate Gaussian halo model, 
a dedicated MCMC analysis of directional data would al- 
low us to constrain the velocity dispersions, resulting in 
a discrimination between various halo models. 



C. WIMP parameters 

As stated above, this MCMC analysis also allows us 
to constrain the parameters of the WIMP by consider- 
ing both the angular and the energy information from 
each recoiling event. Figure [3] (flrst 2 columns) presents 
marginalized distributions and 2D correlation plots con- 
cerning the WIMP parameters (m^, logiQ(CT„)). First, 
we can notice that this analysis method allows us to get 
satisfactory results, i.e. constraints which are consistent 
with the input values and with a rather small dispersion: 

= 51.8tl9,i GeV/c2 (68% CL), 
logio(a„) = -S.Oltg.^s (g8% CL) 

Moreover, as the velocity dispersions are set as free 
parameters, induced bias due to wrong halo model 
assumptions is avoided as long as the input halo model 
is consistent with our ansatz. We refer the reader to 
[Tsj for a detailed discussion about the effect of halo 
model uncertainties on allowed regions. In fact, the 
combined use of angular and energy information allows 
us to remove degeneracies amongst the eight parameters 
and hence to obviate bias in the determination of the 
WIMP properties. 

We observe the usual strong correlation 
/9[m^, logio(o'n)] ~ 1 (see Fig. \^ between and 
log2^Q(fT„) which is inherent in the very definition of 
the event rate, as it scales basically with a„/m^ for 
low mass target. We also found a small and positive 
correlation between the WIMP mass (inversely propor- 
tional to /is) and the background rate (proportional to 
Pb) such as p[m-^,Rb] ~ 0.25. Indeed, as mentioned 
before, this correlation is straightforwardly due to the 
relationship between /is and pb where the total number 
of recorded events follows a Poisson distribution of mean 
Ps + Pb ~ Ncvcnt- Finally, we found no correlation 
between log the background rate Rb. 

As a conclusion, directional detection provides a 
unique opportunity to constrain, with a single experi- 
ment, the WIMP mass and the WIMP-nucleon cross sec- 
tion within the framework of a high-dimensional multi- 
variate analysis. This is of grea t interest in the context 
of phenomenological efforts [Tl - flTl [l9l 1451 - 147} trying to 
constrain the WIMP parameters (m^, ct„) with upcoming 
dark matter experiments, with either an indirect, direct, 
or directional strategy. In this work, we have gone one 
step further in constraining the local WIMP velocity dis- 
tribution with a single directional detection experiment. 
It is, of course, possible to include external data as nui- 
sance parameters, e.g. measurement of the local dark 
matter density po [48, 49,], the local circular velocity 
vq and the escape velocity (taken as infinity in this 
study). However, it seems premature at the level of 
a methodological study aimed at showing how to han- 
dle directional detection data. For instance, the lo- 
cal WIMP density is usually quoted within the range 
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po ~ 0.2 — 0.8 GeV.c~^.cm~'^ and we used the so-called 
"standard" value 0.3 GeV.c~^.cni~^ for the sake of com- 
parison with various direct detector results [il] . We note 
that recently, a value of the local dark matter density, 
Po = 0.43 ± 0.11 ± 0.10 GcV.c-2.cm-3^ has been derived 
within the framework of a Galaxy-model-independent 
method [ioj . All constraints on the WIMP-nucleon cross 
section can be relaxed into constraints on po x a. 



D. Background estimation 

The background rate estimation is also a key point 
of this analysis strategy, not for the value itself but for 
the fact that a wrong background estimation may induce 
bias for other parameters. Indeed, as upcoming data will 
necessarily be contaminated by some background events, 
it is important to be able to manage them. As shown 
in Fig. [3] (last row), it is correctly estimated from the 
MCMC : Rt = 10.97 ± 1.2 kg-^year^^ (68% CL), with 
tiny correlations with other parameters already discussed 
in the previous sections. Then, the fact that the back- 
ground rate is left as a free parameter and reconstructed 
with the MCMC method allows us to avoid bias in the 
estimation of the other parameters. Qualitatively, the 
background rate is mainly constrained by the angular 
part of the spectrum, more precisely in the hemisphere 
opposite to the Cygnus constellation, where few WIMP 
events are expected. In fact, the quality of the estimation 
of the WIMP and halo parameter is directly related to the 
estimation of the background rate. In this example, we 
have shown that dark matter parameter estimation (main 
direction, WIMP, and dark matter halo properties) is not 
affected by a rather large background fraction (^ 30%). 
Hence, directional detection can accommodate to a size- 
able background contamination (posterior to data selec- 
tion), suggesting the idea that light shielding might be 
sufficient, thus allowing us to reduce muon-induced neu- 
tron background (50| . 

As stated above, for this example a fiat background en- 
ergy spectrum has been considered, which is indeed an 
optimistic case. In Sec. |Vll we study the effect of con- 
sidering an energy distributions for background events 
which is similar to the one for WIMP events. 



VI. RESULTS FOR VARIOUS INPUT MODELS 

The constraints on the different parameters obviously 
depend on the input model, characterized by the 
WIMP and dark matter halo properties as well as the 
background energy spectrum. Indeed, the directional 
WIMP event rate crucially depends on the dark matter 
parameters, both from particle physics and Galactic halo 
physics, and degeneracies may arise depending on their 
input values. In the following, we explore various input 
models in order to evaluate their impact on the different 
constraints which could be obtained with a single 
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FIG. 6: 95% contour level in the {£q, 60) plan for three input 
models: Isotropic halo model 4- exponential background (solid 
line). Isotropic halo model + flat background (long dashed 
line) and anisotropic halo model + flat background (dotted 
line). 



directional detection experiment, as the one proposed 
by the MIMAC collaboration, using our MCMC analysis. 

The first point worth emphasizing is the fact that in all 
cases presented hereafter, the recovered main recoil direc- 
tion is always pointing towards Cygnus, within at most 
~ 4° at 95% CL (see Fig.E]). This is relatively straight- 
forward, given the fact that this directional signature is 
uncorrelated with the other parameters of the MCMC 
analysis, as emphasized in Sec. IV Al Indeed, it has been 
shown in [1] that this directional signature only depends 
on the background contamination, which is taken equal 
to 10 evts/kg/year in every following cases. This out- 
lines the robustness of the choice of this parameter as a 
relevant observable to prove that a positive detection of 
dark matter has been reached by a directional detector. 
As outlined in Q, this would allow directional detection 
to provide evidence in favor of a detection of Galactic 
dark matter even at low exposure and even with a size- 
able background contamination. In this study, we have 
checked that this conclusion holds true even in the case 
of non standard dark matter halo model. 



A. Varying the input WIMP mass 

As highlighted by several previous studies p^ - [l^ [l9| , 
the WIMP mass plays a key role in the shape of the 
allowed regions. We have simulated three different sets 
of directional data corresponding to an input WIMP 
mass of = 20,50,100 GeV/c^ with a constant 

WIMP-nucleon cross section an — 10^'^ pb, considering 
a MIMAC-like directional detector (Sec. ImBl) and the 
standard isotropic halo model. The results from the 
three MCMC runs are illustrated in Fig. [71 We present 
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FIG. 7: Left panel : 68% and 95% contour level in the {mx,crn) plane, for the isotropic input model and for a WIMP mass 
equal to 20, 50 and 100 GeV/c^. Right panel : posterior PDF distribution of the /3 parameter for the same models. 




FIG. 8: Left panel : 68% and 95% contour level in the (m^,(T„) plane, for a 50 GeV/c^ WIMP and for two input models : 
isotropic (/? — 0) and triaxial (/3 — 0.4). Right panel : posterior PDF distribution of the /3 parameter for the same models. 
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FIG. 9: Left panel : 95% contour level in the {m^,an) plane, for a 50 GeV/c^ WIMP and for three input background model : 
no background, flat spectrum and exponential spectrum. Right panel : posterior PDF distribution of the /3 parameter for the 
same models. 
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TABLE III: Values of the fi parameter from the marginalized 
distribution for various input models. We quote mean value 
of the PDF distribution and 68 % CL error bars. 



for the three WIMP masses, in the left panel, the 68% 
and 95% CL contours in the (m^, log]^o(f''n)) plane, and 
in the right panel, the posterior PDF P(/?|Z?) of the 
anisotropy velocity parameter /?. 

The WIMP properties (mp,,, logjQ(a-„)) are consistently 
constrained according to the input values with no a 
'priori knowledge of the halo properties, as the velocity 
dispersions are set as free parameters of the analysis. It 
can be deduced from Fig. [7] that this analysis is working 
for any input WIMP mass even if the constraints 
strongly depend on the input value. Indeed, as it can 
be seen in Fig. [3 the constraints on (m^, log;iQ(cr„)) 
are very tight below 50 GeV/c^ and become wider for 
increasing WIMP mass. 

In fact, for the 100 GeV/c^ input WIMP mass, only 
a lower limit should be deduced as > 30 GeV/c^ 
(68% C.L.). Indeed, the 68% and 95% CL. contours 
correspond to the case where a flat prior on £ [5, lO'^] 
GeV/c^ is considered. These weaker constraints in the 
case of a heavy WIMP are due to the fact that the signal 
characteristics, i.e the slope of the energy distribution 
and the width of the angular distribution, evolve slowly 
with the WIMP mass once > 100 GeV/c^ for a 
fluorine target and a recoil energy in the range [5,50] 
keV, as shown in Q. 

As a consequence of this weaker constraint at heavy 
WIMP masses, the constraints on the halo properties are 
also getting weaker, albeit with smaller effect. Indeed, 
as shown in the right panel of Fig. [71 the constraint on 
the anisotropy velocity parameter /? is stronger (smaller 
error bars) for an input WIMP mass of 20 GeV/c^ 
than for a 100 GeV/c^ one. However, as highlighted in 
Table IIIIl the constraint on the /3 parameter remains 
competitive and for a 30 kg. year exposure with a 
MIMAC like directional detector, this MCMC would 
allow us to get, in this strong evidence in favor 

of an isotropic dark matter halo. 



B. Effect of an anisotropic input halo model 

In this section, we vary the input halo model to 
evaluate the evolution of the constraints associated with 
the different dark matter properties (m^^, (T„, /?). Indeed, 
as the velocity dispersions are set as free parameters, 
induced bias due to wrong model assumption should be 
avoided. This is for instance the effect observed in [l3j . 
with a systematic downward shift of the estimated cross 
section, when assuming a standard isotropic velocity 
distribution fitting model whereas the input model is a 
triaxial one. 

In the following, we investigate the effect of an extremely 
triaxial input halo model with (3 = 0.4 (see sec|ll]) on 
the estimation of the dark matter parameters. 

The results from the MCMC run on a simulated 
dataset corresponding to a WIMP mass of 50 GeV/c^ 
with the latter anisotropic halo model are presented 
in Fig. [HI As for the previous section, in the left 
panel is presented the constraint at 68% and 95% on 
the (m^, log]^Q(fT„)) plane, while in the right panel is 
given the deduced posterior PDF of the /? parameter. 
For convenience and comparison, the results from the 
benchmark input model (isothermal sphere with a 50 
GeV/c^ WIMP) are recalled. 

From the left panel of Fig. [SJ we can conclude that 
the two halo models give similar constraints which are 
both consistent with the input values. In fact, and 
as foreseen, the fact that the velocity dispersions are 
set as free parameters in the MCMC analysis allows 
us to avoid induced bias due to wrong model assumption. 

From the right panel of Fig. |S] we can deduce that the 
(3 parameter is well constrained: /? — 0.38lg ?, as in the 
isotropic case. In fact, the constraint is even stronger 
in the anisotropic case than in the isotropic one. This 
comes straightforwardly from the decrease of the degen- 
eracy between the three velocity dispersions with increas- 
ing departure from isotropy. 

As a conclusion of this study, it should be highlighted 
that the combination of information from the angular 
and energy distributions leads to robust allowed regions 
in the (m^, log]^Q((7„)) plane, since the halo model is 
also being constrained with the MCMC analysis from 
the same dataset of a single directional detection exper- 
iment. Moreover, the velocity anisotropy parameter /3, 
i.e. the three velocity dispersions, could be sufficiently 
constrained to discriminate between different halo models 
with future directional detectors such as the one proposed 
by the MIMAC collaboration [1]. 

C. Varying the input background spectrum 

The background energy spectrum is a key issue 
for both directional detection and direct detection 
(direction-insensitive experiments). When setting ex- 
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elusion limits with directional detection, the difficulty 
can be avoided by considering only the angular part 
of the directional event rate, thus allovifing us to set 
robust and conservative limits Q- But as far as the 
whole directional event rate is used, the question of the 
background energy spectrum must be carefully treated. 
In fact, as the background energy spectrum is unknown 
it must be guessed to be included in a likelihood-type 
analysis. Then, a wrong assumption on the background 
shape leads to an incorrect estimation of the background 
rate, resulting in a wrong estimation of the dark matter 
properties. 

Motivated by simulations of neutron background in dark 
matter detectors, e.g. in low pressure TPC [5l| . two 
different background energy distributions are usually 
considered |12h17| : flat and exponentially decreas- 
ing with increasing recoil energy. The exponential 
one corresponds to the most pessimistic case as it is 
chosen, in our case, to be exactly the same as the 
WIMP-induced energy distribution. That is to say 
an exponential distribution with a slope of ^ 17 keV 
in the case where the WIMP mass is 50 GeV/c^ and 
considering an isotropic halo model. As outlined in [l^ . 
it is not possible to disentangle a WIMP signal from the 
background, with a single direct detector, if the shape of 
the background and WIMP-induced energy distributions 
are similar. In principle, this will not be the case for 
directional detection as the angular distribution of the 
background is isotropic, then remaining different from 
the WIMP-induced one. 

Figure[H]presents the constraints in the (m^, logig(cr„)) 
plane and on the /3 parameter for three input background 
energy distributions: no background (black solid line), 
flat (red dashed line) and exponential (blue dotted line). 
The other input parameters are those used in Sec. |V] 
for the benchmark input model (an isotropic halo, a 
50 GeV/c^ WIMP mass and a 10"^ WIMP-nucleon 
axial cross section). 

First, the comparison between the case with no back- 
ground and the case with a flat background energy 
distribution highlights the fact that even with a large 
background contamination (~ 30%), the results are quite 
similar, particularly in the determination of the WIMP 
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properties, due to the fact that the disentanglement 
between WIMP and background events is done with 
both energy and directional arguments. Then, in both 
cases, as the background rate is correctly estimated by 
the MCMC analysis (see Table lIVp , systematic bias in 
the estimation of the dark matter properties is avoided. 
However, from Fig. [S] and Table Hill it should be noticed 
that, even if the /3 parameter is consistently constrained 
according to the input value, the presence of a sizeable 
background leads to a wider constraint (about two times 
larger). 

In the case of an exponential input background energy 
distribution, the result is basically unchanged, although 
constraints are weaker. This is due to the fact that the 
background event rate parameter is less constrained (see 
Table |IV|) resulting in broader marginalized distributions 
of other parameters. Indeed, the estimation of the Ri, 
parameter is done solely with the angular part of the 
spectrum, as the energy distributions are exactly the 
same for both kinds of events. Nevertheless, even in such 
a pessimistic case, the WIMP properties (m^, logjQ(f7„)) 
and the dark matter halo properties encoded in the 
/3 parameter can still be estimated with upcoming 
directional detectors with realistic exposures thanks to 
the use of this MCMC analysis. 

From this study, it should be concluded that, the effect 
of background contamination on directional data can be 
handled in the case where the background energy distri- 
bution is correctly estimated. Eventually, we have shown 
that even for a large background contamination and in 
the most pessimistic background model, directional de- 
tection combined with this MCMC analysis should allow 
us to assess consistent and interesting constraints on the 
dark matter properties with a single experiment. 



VII. CONCLUSION 

We have shown that identification of dark matter 
might be achieved with a 10 kg CF4 directional detector, 
allowing 3D track reconstruction with sense recognition 
down to 5 keV and operated during three years. To fully 
exploit upcoming data, we propose a new high dimen- 
sional multivariate analysis method based on a Markov 
chain Monte Carlo analysis of recoil events, allowing to 
constrain, in a single directional experiment, the WIMP 
parameters, both from particle physics (mass and cross 
section) and Galactic halo (velocity dispersion along the 
three axis) and within the framework of a given ansatz. 

Indeed, the combination of information from the angu- 
lar and energy distributions leads to robust allowed re- 
gions in the (m^, log]^Q(cr„)) plane, since the halo model 
is also being constrained with the MCMC analysis from 
the same dataset of a single directional detection exper- 
iment. Moreover, the velocity anisotropy parameter /3, 
related to the three velocity dispersions, could be suffi- 
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ciently constrained to discriminate between various halo methods, 
models with future directional detectors such as the one 
proposed by the MIMAC collaboration Q. 
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